#include<stdio.h>
#include<math.h>
#include<time.h>

#define SAMPLE_SIZE 2000
#define PI 3.14159265358979f
void writeDataFile(float[], int);

int MT[624];
int index = 0;

void init_gen(int seed){
	MT[0] = seed;
	for(int i = 1; i < 624; i++)
	{
		MT[i] = 1812433253 * (MT[i-1] ^ (MT[i-1]>>30)) + i;
		MT[i] &= 0xffffffff;
	}
}

 void generate_number(){
	 int y;
	for (int i = 0; i < 624; i++)
	{
		y = (MT[i] & 0x80000000) + (MT[(i+1) % 624] & 0x7fffffff);
		MT[i] = MT[(i + 397) % 624] ^ y >> 1;

		if((y % 2) != 0)
			MT[i] ^= 0x9908b0df;
	} 
}

float extract_number(){

	if(index == 0)
		generate_number();

	int y = MT[index];
	
	y ^= y >> 11;
	y ^= (y << 7) & 0x9d2c5680;
	y ^= (y << 15) & 0xefc60000;
	y ^= y >> 18;
   
  index = (index + 1) % 624;

	//return y;
	float x = (float)y/0xffffffff;
	return x;
	
}

void boxMuller(float *n1, float *n2){

	float r = sqrt(-2 * log(*n1));
	float phi = 2 * PI * (*n2);

	*n1 = r * cos(phi);
	*n2 = r * sin(phi);
}

int  main(){

	const int sampleSize = SAMPLE_SIZE;
  
  time_t now = time(0);

  // only time(x) works on both platforms.
//	SYSTEMTIME sTime;
//	GetSystemTime(&sTime);
//	long seed = now * 1000 + sTime.wMilliseconds;
   
	init_gen(now);
	float num[sampleSize];
	
	printf("\nThe Random numbers are:\n");
	for(int i = 0; i < SAMPLE_SIZE; i++)
	{
		num[i] = extract_number();
		printf("%f, ", num[i]);
	}

	printf("\n\nAfter box muller transformation:\n");
	for(int i = 0; i < sampleSize / 2; i++)
	{
		boxMuller(&num[i], &num[i + (sampleSize / 2)]);
		printf("%f, %f, ", num[i], num[i + (sampleSize / 2)]);
	}

	//writeDataFile(num, sampleSize);
   
	getchar();
	return 0;
}

void writeDataFile(float data[], int length){
	FILE *f1;
	
	f1=fopen("data.txt","w");
	fprintf(f1, "#Data\n");

	for(int i = 0; i < length; i++){
		printf("%i %f\n", (i+1), data[i]);
		fprintf(f1, "%i, %f\n", (i+1), data[i]);
	}

	fclose(f1);
}
